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The ST2 interaction potential has been used in a large number of simulation studies to explore the 
possibility of a liquid-liquid phase transition (LLPT) in supercooled water. Using umbrella sampling 
Monte Carlo simulations of ST2 water, we evaluate the free energy of formation of small ice nuclei 
in the supercooled liquid in the vicinity of the Widom line, the region above the critical temperature 
of the LLPT where a number of thermodynamic anomalies occur. Our results show that in this 
region there is a substantial free energy cost for the formation of small ice nuclei, demonstrating that 
the thermodynamic anomalies associated with the Widom line in ST2 water occur in a well-defined 
metastable liquid phase. On passing through the Widom line, we identify changes in the free energy 
to form small ice nuclei that illustrate how the thermodynamic anomalies associated with the LLPT 
may influence the ice nucleation process. 


I. INTRODUCTION 

The proposal that a liquid-liquid phase transition 
(LLPT) occurs in supercooled water has stimulated a sig¬ 
nificant quantity of research and debate [Tj-filil . The prin¬ 
cipal virtue of this proposal is that it provides a common 
origin for two of water’s most unusual properties. First, 
the atypical and increasingly divergent behaviour of ther¬ 
modynamic and dynamical properties of liquid water as 
temperature decreases through the freezing temperature 
can be understood in terms of an approach to the critical 
point of the LLPT, and its associated anomalies. Second, 
the line of first-order phase transitions that terminates at 
the critical point of the LLPT can account for the occur¬ 
rence of two distinct glass forms of water, the so-called 
low-density amorphous and high-density amorphous ices. 

The LLPT proposal predicts the occurrence of an in¬ 
terrelated set of thermodynamic features in supercooled 
water. These features are illustrated in the plane of tem¬ 
perature T and pressure P in Fig. [L] which summarizes 
the behaviour reported in previous molecular dynamics 
(MD) simulations of the ST2 model of water (si. [l3|. The 
coexistence line separating the distinct low-density liquid 
(LDL) and high-density liquid (HDL) phases has nega¬ 
tive Clapeyron slope and terminates at a critical point at 
T = T c . Emanating from the critical point for T < T c , 
and bracketing the coexistence line, are two metasta¬ 
bility limits (spinodals), one for each of the LDL and 
HDL phases. For T > T c , a number of thermodynamic 
response functions (e.g. the isothermal compressibility 
Kt) exhibit maxima which diverge as T —> T c . The 
loci of these maxima in the (T, P) plane all lie near the 
so-called Widom line Q, the locus of maxima of the cor- 
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relation length, and converge at the critical point. The 
region near the Widom line therefore contains the high- 
temperature harbingers of the critical point and the line 
of first-order phase transitions that occur for T < T c . 
Whereas the HDL-LDL transition is discontinuous at 
the line of first-order transitions, the liquid continuously 
crosses over from an HDL-like liquid to an LDL-like liq¬ 
uid as T decreases through the region of the Widom line. 

The main source of challenge and concern related to 
the LLPT hypothesis is ice formation. Experimentally, 
rapid ice formation has so far thwarted the direct mea¬ 
surement of the properties of bulk supercooled liquid wa¬ 
ter in the region of T and P where the LLPT is predicted 
to occur [Si]. That is, neither the first-order LLPT, the 
critical point, nor the response function maxima associ¬ 
ated with the Widom line have been directly observed in 
experiments on bulk samples of supercooled liquid water. 
However, a number of indirect experimental results are 
consistent with the LLPT proposal. Examples include 
experimental studies of metastable ice mel ting lines Q, 
the behaviour of the amorphous ices [lS], l23| . and the 
thermodynamic and transport properties of confined liq¬ 
uid water [l6| . Recent experiments continue to push 
deeper into the supercooled liquid regime, offering new 
insights [HI] . 

To date, it is only in computer simulations of water¬ 
like models that a LLPT has been unambiguously ob¬ 
served [24|, [2^| . The most extensively studied model in 
this regard is the ST2 interaction potential [13 ■ In ST2 
simulations, most previous work shows that it is possible 
to evaluate the equilibrium properties of the deeply su¬ 
percooled liquid on a time scale that is short compared 
to the crystal nucleation time, thus expo sing the physics 
of the LLPT to direct observation [25|, |29|. However, 
Limmer and Chandler have disputed this interpretation 
of the results for ST2 water, and have specifically pro¬ 
posed that the behaviour previously interpreted as ev- 
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FIG. 1: Phase behaviour of the ST2 model predicted from pre¬ 
vious MD simulations S M- Shown are the estimated loca¬ 
tions of the liquid-liquid critical point (circle with error bars); 
the HDL spinodal (down-triangles) and the LDL spinodal (up- 
triangles); the line of Kt maxima (thick solid line) and min¬ 
ima (thick dashed line); the line of density maxima (thin solid 
line) and minima (thin dashed line); and the metastability 
limit of the liquid at negative pressure (diamonds). The dot- 
dashed line is a line having the estimated slope of the liquid- 
liquid coexistence curve at the critical point. The circles la¬ 
belled A and B locate the state points at which we carry out 
the umbrella sampling simulations described in the text. The 
enclosed polygon delimits the region of states within which 
we estimate the free energy of formation of small ice nuclei. 
The contours within this polygon are lines along which G 15 is 
constant. The contours are lkT apart, and range from 22 kT 
in the lower left corner to 46 kT in the upper right. 


idence of a LLPT is in fact a manifestation of crystal 
formation [HI, HH H3] • Subsequent studies by others us¬ 
ing very similar methods do not reproduce the simulation 
results of Limmer and Chandler, and confirm the occur¬ 
rence of a LLPT in ST2 water [ii| Hfj- HH, [25| . The origin 
of this discrepancy remains a subject of debate [13, . 

In light of the above results, the nature of crystal for¬ 
mation in a liquid exhibiting a LLPT, and in ST2 water 
in particular, merits further investigation. The relation¬ 
ship between ice formation and the potential for a LLPT 
varies widely in water-like computer models. For exam¬ 
ple, in the mW model of water, recent work shows that 
the supercooled liquid becomes unstable to ice formation 
before the conditions for a LLPT can be realized [H, [lH . 
At the other extreme, simulations of a water-like patchy 
colloidal liquid, and a model of DNA tetramers, demon¬ 
strate that a LLPT can occur in the stable liquid region, 
where crystal formation is impossible [9, (2J, |30j, [3lj . The 
ST2 model provides a case where the physics of a LLPT 
and of ice nucleation occur simultaneously, allowing the 
influence of the LLPT on ice formation to be explored 
directly. The coincident appearance of both phenom¬ 


ena can be understood in terms of the effective rigidity 
of hydrogen bonds when viewed against a broader spec¬ 
trum of bond flexibility in tetrahedral network-forming 

liquids [23.130,. 36]. 

In the most recent work on ST2 water, the free energy 
surface of the system near the critical point of the LLPT 
was examined as a function of the density p and the bulk 
bond-order parameter Qe- using the method of umbrella 
sampling [lj fl9l-l22l. (33 . Q e is a bulk measure of the 
overall crystallinity of a system developed by Steinhardt 
et ah; see e.g. the definition of Qe given in Ref. a As 
such Q 6 is a valid order parameter for distinguishing be¬ 
tween the disordered liquid and ordered crystal phases, 
and provides a reaction coordinate to quantify the varia¬ 
tion of the free energy associated with the crystal-liquid 
phase transition. 

We focus here on the related question of how single ice 
nuclei form and grow from the supercooled liquid. The 
free energy G(n) to form an ice-like cluster of n molecules 
from the liquid can be evaluated from computer simula¬ 
tions using, 

{3G(n) = - hi^^, (1) 

where /? = 1/kT , k is Boltzmann’s constant, N is the 
number of molecules in the system, and J\f(n) is the equi¬ 
librium average of the number of ice-like clusters of size 
n (37j. So defined, G(n) can be used to determine the 
nucleation free energy barrier appropriate for estimat¬ 
ing the nucleation rate, e.g. via classical nucleation the¬ 
ory p38l ]. When using umbrella sampling methods to com¬ 
pute W(n), recent work has shown that an appropriate 
order parameter to define a path from the liquid to the 
crystal phase is n max , the size of the largest crystalline 
cluster in the system [39l - [4ll |. In particular, Ref. [4fj 
discusses why the order parameter rc max is preferable to 
order parameters such as Qe for evaluating the free en¬ 
ergy of formation of crystal nuclei when the aim is to 
estimate the nucleation barrier and rate. 

In order to fully elucidate the influence of a LLPT 
on crystal nucleation, it would be valuable to determine 
G{n) for ST2 water near the critical point of the LLPT, 
and for a sufficient range of n to find the nucleation free 
energy barrier and the size of the critical nucleus. How¬ 
ever, in order to achieve this goal, a number of concerns 
present themselves. For example, G(n) has not been pre¬ 
viously reported for the ST2 model, and has been found 
for only a few other water-like models when using n max 
as the order parameter [ 3 ^, [3l| . An initial assessment of 
the methods for finding G(n) for ST2 water is therefore 
warranted. In addition, relaxation times near the critical 
point of the LLPT are large, relative to easily accessed 
computational time scales, and it is therefore appropriate 
to work first in a regime where equilibration time scales 
are readily attained. 

Consequently, in the present work, we evaluate G(n) 
in ST2 water for small ice nuclei in the range 0 < n < 20, 
and focus on T and P in the vicinity of the Widorn 
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FIG. 2: Representative configurations of ice Ic, ice Hi, and 
liquid water obtained from MC simulations at T = 280 K 
and P = 40 MPa. The simulations of ice Ic and the liquid are 
conducted in a cubic simulation cell with N = 216 molecules. 
The ice Ih simulation is carried out in a rectangular cell with 
N = 432 molecules. 


line. These choices allow us to test the methods for 
estimating G(n ) for ST2 water using a relatively small 
system (N = 216 molecules) and in a range of T and 
P where equilibrium is easily established. Further, for 
finding G(n), we explore the use of umbrella sampling 
with respect to the two order parameters p and n max . 
As described below, umbrella sampling with respect to p 
allows us to reweight our results with respect to P over a 
range near the chosen P of our simulations. We also con¬ 
duct separate series of runs at two values of T, and are 
able to interpolate results in between by reweighting with 
respect to T. The result is an estimate of G(n) inside a 
continuous region of states in the (T, P) plane straddling 
the Widorn line. Thus, although the present study does 
not explore the immediate vicinity of the critical point 
of the LLPT, our results allow us to draw conclusions 
about the stability of the metastable liquid state near 
the Widorn line, and to discern how the thermodynamic 
anomalies that extend outward from the LLPT may in¬ 
fluence ice nucleation. 


II. MODEL AND ORDER PARAMETERS 



FIG. 3: Normalized probability distributions P(dj) repre¬ 
senting the frequency of occurrence of dj values for the equi¬ 
librium liquid, ice Ic, and ice Ih phases at T = 280 K and 
P = 40 MPa. 



We study the ST2 model of water proposed by Still- 
inger and Rahman p33j ]. The details of our implementa¬ 
tion of the ST2 interaction potential are summarized in 
Ref. |2C}, and are the same as those used in a number of 
previous studies 012,! EMI- In particular, we note 
that in our implementation of the ST2 model, we use 
the reaction field method to approximate the long-range 
contribution of the electrostatic interactions [42). 

To characterize the crystallinity of the system, we fol¬ 
low the approach of Refs. [H, [43[. We use Steinhardt 
bond order parameters [44| based on spherical harmon¬ 
ics Yim of order 1 = 3. For each particle i, and for integers 
to such that —l < m < l , we define the components of a 
complex vector as, 


FIG. 4: Normalized probability distributions P(ribonds) for 
the number of crystalline bonds nbonda found for molecules in 
the equilibrium liquid, ice Ic, and ice Ih phases at T = 280 K 
and P = 40 MPa. 

where the sum is over the four nearest neighbours of 
molecule i, and fik is a unit vector pointing from the 
O atom of molecule i to that of molecule k. The dot 
product, 

i 

Cij= ( 3 ) 

m=—l 

where, 
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and q* m {i) is its complex conjugate, quantifies the orien¬ 
tational correlation between the nearest-neighbour shells 
of molecules i and j. 

We carry out initial simulations of ST2 water in the liq¬ 
uid phase, and in the ice Ic and Ih crystalline phases, all 
at the same state point at T = 280 K and P = 40 MPa; 
see Fig. [2] Fig. [3] shows the normalized probability distri¬ 
butions of Cij for each phase. The distribution for the liq¬ 
uid phase is broad and spread out over much of the range 
of c^. The ice Ic distribution is peaked sharply around 
Cij = — 1, reflecting the fact that each molecule partic¬ 
ipates in four “eclipsed” hydrogen bonds in the ice Ic 
structure. The ice Ih distribution exhibits two peaks at 
approximately Cij = —1 and c,; ? = —0.1, indicating the 
presence of three eclipsed, and one “staggered” hydrogen 
bond for each molecule in this crystal. 

Following Refs. [H, |42j, we define that a crystalline 
bond exists between molecule i and one of its four nearest 
neighbours j when Cy < —0.87. A crystalline molecule is 
defined as one that has three or more crystalline bonds. 
Fig. H] demonstrates that this definition for a crystalline 
molecule is unlikely to be met by molecules in the liquid 
phase. Conversely, this definition is unlikely to misiden- 
tify molecules in either crystal phase as liquid-like. 

As discussed in Refs. this definition permits 

the formation of both the ice Ic and Ih structures, since 
both have at least three eclipsed bonds. As a conse¬ 
quence, the umbrella sampling procedure that we will 
apply to the liquid state will not specifically select for 
one of these crystal phases to the exclusion of the other. 
At the same time, the possibility exists that this approach 
favors the ice Ic structure over that of ice Ih, since only 
three of the four eclipsed bonds of an ice Ic molecule 
need to be detected, while all three of the eclipsed bonds 
of an ice Ih molecule must be detected. In the study of 
a patchy colloid model Jijj, the influence of this bias to 
form ice Ic was indeed observed, although the height of 
the nucleation barrier was not affected. 

Having identified all the crystalline molecules in the 
system, we can define crystalline clusters. If a crystalline 
molecule j is one of the four nearest neighbours of a crys¬ 
talline molecule i. then i and j are assigned to the same 
crystalline cluster. The size n of a crystalline cluster is 
the total number of crystalline molecules contiguously 
connected in this way. We then find n max , the size of the 
largest crystalline cluster in the system. We use n max as 
the order parameter to describe the progression of the 
system from the liquid to the crystalline phase. 

III. UMBRELLA SAMPLING 

We carry out umbrella sampling MC simulations in the 
constant- (TV, P, T) ensemble. Our approach is similar to 
that used in Refs, [lj, []j| Ho[, except that we use n m ax 
rather than Qq as a structural order parameter. To im¬ 
plement umbrella sampling, we add a biasing potential, 

U B = hip-p*) 2 + fc 2 (n max - n* max ) 2 , (5) 



FIG. 5: Time dependence of the autocorrelation functions 
C x {t ) for x = n ma x, U, and p, as obtained from the series A 
run with p* = 0.80 g/cm 3 and n^ ax = 19. This the slowest 
relaxing run of series A. 


to the system potential energy U. The effect of Ub is 
to constrain a given simulation to sample configurations 
in the vicinity of chosen values of the order parameters 
p = p* and n max = rt* nax . In all our simulations, we fix 
N = 216, h = 1000 kT (cm 3 /g) 2 , and k 2 = 0.5 kT. 

Trial configurations for each Monte Carlo step (MCS) 
are generated as follows: First, we carry out an unbiased 
(i.e. Ub = 0) constant-(W P, T) MC “sweep”, in which 
each we attempt (on average) N — 1 rototranslational 
moves, and one change of the system volume. The maxi¬ 
mum size of the attempted rototranslational and volume 
changes are chosen to give MC acceptance ratios in the 
range 25-40%. Next, the change in the biasing potential 
Ub is evaluated for the trial configuration resulting from 
the unbiased MC sweep, relative to the system configu¬ 
ration at the beginning of the sweep, to determine the 
acceptance or rejection of the trial configuration. This 
completes one MCS, and the procedure is then repeated. 

We next identify the (T, P) state points near the 
Widom line at which to conduct our runs. As a proxy for 
the Widom line, we use the line of compressibility max¬ 
ima (LCM) shown in Fig. [TJ as determined from previous 
MD simulations of the ST2 model [§|. We select two state 
points lying on the LCM at relatively high T, in order 
to study conditions at which the relaxation time of the 
liquid will be relatively short. In the following, “series 
A” denotes the set of runs conducted at T = 300 K and 
P = — 25 MPa, and “series B” denotes the set of runs 
conducted at T = 280 K and P = 40 MPa. The state 
points in the (T, P) plane corresponding to series A and 
B are identified in Fig. [0 For both series, we carry out 
210 independently initialized umbrella sampling runs at 
equally spaced values of p* from 0.80 to 1.20 g/cm 3 sep¬ 
arated by 0.02 g/cm 3 ; and equally spaced values of n^ax 
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FIG. 6: Time dependence of the autocorrelation function Cu(t ) for series A and B runs having n^ax = 1 and 19. In each panel, 
the curves are for p* = 0.8, 0.9, 1.0, 1.1, and 1.2 g/cm 3 , from most slowly decaying to most rapidly decaying. 


from 1 to 19 separated by 2. 

To initialize our umbrella sampling runs, we first con¬ 
duct 210 unbiased MC runs at state point A starting 
from independent configurations harvested from a high- 
temperature run conducted at T = 600 K and P = 
25 MPa. We separately equilibrate all of these unbiased 
runs at state point A. Then the umbrella sampling po¬ 
tential Ub is imposed, and the time evolution continued 
in order to bring the biased runs into equilibrium. In all 
cases we find that a crystalline cluster of size n max close 
to n \ nax develops during each run when initiated from 
the liquid state. After all of the runs in series A attain 
equilibrium, we use a configuration from each to initi¬ 
ate the corresponding umbrella sampling run for series B 
by changing the imposed T and P. These series B runs 
are then separately evolved to equilibrium. From each 
run, we record a time series of instantaneous values for 
p, nmax, U, and for Af(n) for n = 0 to 20. 

Our runs are carried out for between 4.5 x 10® and 
6.3 x 10 8 MCS. Using the second half of each run, we com¬ 
pute the autocorrelation function C x (t) from the time se¬ 
ries for each choice of x = p, n max , or U. As illustrated 
in Fig. 0 we find in almost all cases that U is the most 
slowly relaxing of these three quantities. As shown in 
Fig. [6l in all cases Cu(t) decays to zero on a time scale 
that is short compared to the lengths of our runs, con¬ 
firming that our runs are well equilibrated. To charac¬ 
terize the relaxation time of our umbrella sampling runs, 
we choose the time r such that =0.2. As shown 

in Fig. [Tj we find values of r ranging between 9 x 10 8 and 
3 x 10 6 MCS. To account for equilibration, we then dis¬ 
card the results for t < T e of each run, where r e is chosen 
to be at least 20r. The resulting length r run of each pro¬ 
duction run that is used in our analysis is shown in Fig. [7] 
compared to the corresponding value of r. In terms of t, 
the lengths of our production runs range between 173r 


and 57000 t. 

To estimate observables and their associated error, 
we use the multistage Bennet acceptance ratio (MBAR) 
method j4|| . The MBAR method takes as input the time 
series of the order parameters (p and ?r m ax) and the sys¬ 
tem potential energy U, reweights the statistics obtained 
from each run to remove the effect of the biasing po¬ 
tential, and produces an optimal estimate of the desired 
observable at a specified value of T and P. Similar to his¬ 
togram reweighting methods [46], [47j , the MBAR method 
also facilitates reweighting the configurations sampled 
during our runs with respect to T and/or P, allowing 
the statistics from different state points to be combined 
to produce an estimate of a desired observable at (T, P) 
state points that lie near the conditions at which we carry 
out our simulations. 

For the purpose of estimating average values and their 
error using the MBAR method, we wish to consider only 
those configurations from our runs that are statistically 
independent. We assume that statistically independent 
configurations are separated by r or r run /1000, whichever 
is larger. All other configurations are ignored in our anal¬ 
ysis. For the G(n) data presented in the figures that fol¬ 
low, we find that one standard deviation of error is always 
comparable to or smaller than the symbol size. 

In order to determine the range of T and P within 
which we can reliably use reweighting to estimate the av¬ 
erage values of observables from our data set, we must 
examine the range of p and U sampled by our umbrella 
sampling procedure. That is, we can only evaluate ob¬ 
servables at values of P where our runs have adequately 
sampled microstates having values of p that are typical 
of equilibrium at that value of P. Similarly, we can only 
evaluate observables at values of T where our runs have 
adequately sampled microstates having values of U that 
are typical of equilibrium at that value of T. To pro- 
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FIG. 7: Comparion of r (circles) and r run (squares) as a func¬ 
tion of rimax for all the runs in (a) series A, and (b) series 

B. 


ceed, we first identify in Fig. [5] the mean values of U 
and p for liquid states in the vicinity of state points A 
and B. Fig. [8] shows equilibrium isotherms of U versus 
p obtained from previous MD simulations carried out at 
constant ( N,V,T ), for T = 275 to 305 K We then 
compare these equilibrium isotherms with the values of 
U and p of the microstates sampled in our n* = 1 runs, 
since these are the runs most representative of the equi¬ 
librium liquid, in which crystal nuclei are either absent 
or very small. Fig. [8] shows the mean values of U and p 
sampled in our n* = 1 runs for all values of p *, as well as 
the distribution of sampled values as characterized by the 
standard deviation of the set of U and p values sampled 
during each run. 

The dashed lines in Fig. [5] connect the two highest and 
the two lowest density points from series A and B. We use 
these lines to define the limits of adequate sampling with 



U (kJ/mol) 

FIG. 8: Small data points connected by lines are equilib¬ 
rium isotherms of p versus U as obtained from previous MD 
simulations Jjj. From left to right, the isotherms are for 
T = 275 to 305 K, in steps of 5 K. Also shown are mean 
values of U and p at n* = 1 for each value of p* in series 
A (circles) and series B (squares). The points at the lowest 
p are for p * = 0.80 g/cm 3 , and those at the highest p are 
for p* = 1.20 g/cm 3 . The spread in the values of U and p 
sampled at each value of p* are indicated respectively by the 
horizontal and vertical bars through each point. These bars 
extend one standard deviation above and below the mean, for 
both U and p. The lower (higher) dashed line connects the 
p* = 0.80 g/cm 3 (p* = 1.20 g/cm 3 ) points for series A and 
B. The polygon in Fig. [Q encloses all the ( T , P) values cor¬ 
responding to the data points on the equilibrium isotherms 
that lie between the two dashed lines. 



FIG. 9: Comparison of isotherms of P as a function of p 
evaluated from the present data set using the MBAR method 
(filled symbols), and isotherms obtained from previous MD 
simulations (open symbols) [8]. To facilitate comparison of 
the curves, the data for T = 290 K have been shifted upward 
by 40 MPa, and for T = 275 K by 80 MPa. 































7 



FIG. 10: G{n ) at T = 290 K and P = 0 MPa (open circles), compared to fits to the data of Eq. [6] (dashed line) and Eq. [7] (solid 
line). Also shown are sample configurations in which the molecules belonging to the largest crystalline cluster are rendered 
as solid spheres (red for oxygen and white for hydrogen), and where all other molecules are rendered using a stick model. 
All molecules have been shifted so that the centre of mass of the largest crystalline cluster coincides with the centre of the 
simulation cell. For the configurations at n = 9, 15, and 19, the system has been rotated to highlight the similarity to the ice Ic 
structure shown in Fig. [5] 


respect to p along each of the equilibrium isotherms. This 
in turn defines the range of P within which we evaluate 
observables. That is, for a given T we only consider the 
range of P that corresponds to the range of p found be¬ 
tween these two dashed lines. This criterion determines 
the upper and lower boundaries of the polygon shown in 

Fig.m 

To delimit the range of T examined, our criterion is 
that the equilibrium isotherms in Fig.[8]should fall within 
one standard deviation of the mean values of U sampled 
in our umbrella sampling runs. This criterion is satisfied 
between 275 and 305 K, except for T in the vicinity of 
290 K. However, the 290 K equilibrium isotherm lies be¬ 
tween the series A and B data sets, and is also within two 
standard deviations of both. Consequently, the statistics 
of both series contribute to the evaluation of observables 
near 290 K. As we will see below, we confirm that the 
statistics accumulated in our umbrella sampling runs are 
indeed sufficient to compute observables throughout the 
range of T from 275 to 305 K, thus justifying the high 
and low T boundaries of the polygon shown in Fig. [l] 

To test our ability to evaluate average values from our 
data set over the range of T and P indicated by the 
polygon in Fig. [H we show in Fig. [9] isotherms of P versus 
p, for T and P spanning the width and height of the 
polygon. Fig. [9] demonstrates good agreement between 
the isotherms calculated by using the MBAR method to 
find the average value of p at a given T and P from our 
combined series A and B data sets, and the isotherms 
obtained from previous MD simulations [§j. 


IV. RESULTS 

To find G(n) throughout the (T, P) range within the 
polygon of Fig.[l] we combine the statistics from all runs 
in series A and B. We use the MBAR method to compute 
the average value of Af(n) for n = 0 to 20, on a grid of 
(T, P) points within the polygon, where each point is 
separated by 5 K in T and 20 MPa in P. G(n) at each 
state point is then found from Af(n) using Eq. [L] 

A representative result is given in Fig. 1101 which shows 
G{n) found at T = 290 K and P = 0 MPa, a point very 
close to the LCM. Fig. [TO] shows that there is a large 
free energy cost (up to 30 kT) for the formation of small 
ice nuclei in the vicinity the Kt maximum, demonstrat¬ 
ing that the thermodynamic anomalies associated with 
the Widorn line for ST2 water occur in a well defined 
metastable liquid, for which crystal formation can only 
occur by surmounting a substantial free energy barrier. 

In order to confirm the result presented in Fig. 1101 and 
to test for systematic errors in our implementation of 
the MBAR method (including reweighting), we conduct 
an additional test. As shown in Fig. |TT| we compare 
the G(n) curve from Fig. [TU] (obtained using data from 
both series A and B) with the result obtained from se¬ 
ries A alone, and from series B alone. All three curves 
are within error of one another. This agreement con¬ 
firms that the statistics gathered from series A and B 
have independently converged to equilibrium, and that 
the reweighting of data gathered at either state point A 
or B gives consistent results at T and P intermediate 
between A and B. 
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A only, 300 K, -25 MPa 
A only, reweighted to 290 K, 0 MPa 
B only, 280 K, 40 MPa 
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FIG. 11: Test of the results for G(n) obtained using the 
MBAR method with reweighting to T = 290 K and P = 
0 MPa, using three different approaches: (i) using reweighted 
data from series A only (open circles); (ii) using reweighted 
data from series B only (open squares); and using reweighted 
and combined data from both series A and B (plus signs). 
For comparison, we also show G(n) curves (filled symbols) 
obtained from series A and B without reweighting with re¬ 
spect to T or P. 



FIG. 12: G(n) at T = 290 K and P = 0 MPa (open circles), 
compared to fits to the data of Eq. [6] (dashed line) and Eq. 0 
(solid line). 


Visual inspection of the structure of the larger ice nu¬ 
clei (in the range of n = 15 to 20) suggests that most of 
the nuclei formed in our simulations have the structure 
of ice Ic, as illustrated in Fig. [Tol We have not found any 
examples of ice Ih nuclei. As discussed in Section II, this 
result may be a consequence of the order parameter used 
in our umbrella sampling procedure, or it may indicate 
that for the ST2 model the free energy required to form 
small ice Ih nuclei is higher than for ice Ic, under the con¬ 
ditions studied here. Following the approach of Ref. (3!|, 


0 5 10 15 20 25 

n 

FIG. 13: G(n) at T = 290 K over a range of P. Solid lines 
are fits of Eq. 0 

this question could be resolved in future work by using 
modified order parameters that enforce the formation of 
only ice Ic or ice Ih. 

To test if our results may be used to estimate the height 
of the nucleation barrier G* and the size of the critical 
nucleus n*, we use our data for G(n) to obtain a fit of 
an equation for G(n) the form of which is predicted by 
classical nucleation theory: 

G(n) = —an + bn 2 / 3 , (6) 

where a and b are fit parameters jUj. We also fit a mod¬ 
ified equation that includes an additional term and fit 
parameter c: 

G(n) =—an + bn 2 / 3 + cn 1 / 3 . (7) 

A term of the form of that added in Eq.0has been used to 
account for the curvature of the nucleus-liquid interface 
for simple liquids (49|, as well as to model a core-shell 
morphology of the nucleus in ice nucleation t 4l|. 

The fits of Eq. |6]and [T] obtained for G{n) at T = 290 K 
and P = 0 MPa are shown in Figs. ITOl and Fig. [12] While 
both fitting functions do a reasonable job of representing 
the data in the range of small n , they give widely dif¬ 
ferent estimates for G* and n*. Similar results to those 
shown in Fig. [12] are obtained for state points through¬ 
out the (T, P ) region in which we compute G{n). Given 
the large degree of extrapolation required to estimate G* 
and n* from our data, we do not attempt to analyze these 
quantities further here. 

The variation of G[n) with P at T = 290 K is shown in 
Fig. [T3] For a given n, the value of G{n) increases with 
P. This is in line with expectation. At T = 280 K and 
P = 40 MPa, the density of ice Ic and Ih is 0.862 g/cm 3 , 
whereas the density of the liquid is 0.917 g/cm 3 . At 
higher P, the difference in density between the liquid 
and crystal increases. Hence the local structure of the 
liquid is more tetrahedral and ice-like at lower P, where 
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FIG. 14: Isotherms of G 15 as a function of P. Isotherms are 
5 K apart. The open circles locate the pressure of the Kt 
maximum along each isotherm. 


FIG. 15: Isobars of G 15 as a function of T. T c and P c are 
respectively the estimated T and P of the critical point of the 
LLPT. 


the densities of the liquid and crystal are closer. As P 
increases, the local tetrahedral network of the liquid is 
progressively disrupted, and a more significant local re¬ 
structuring is required to create an ice-like cluster. From 
this point of view, we should expect the free energy cost 
to create ice nuclei to increase with P in the range of 
pressure examined here. 

To characterize the behaviour of G{n) over the en¬ 
tire (T, P) range studied here, we focus on the value 
of G 15 = G(n = 15). Isotherms of G 15 are shown in 
Fig. [id] and isobars in Fig. [15] The two dimensional 
function Gis(T, P) is represented in Fig. |T] as a con¬ 
tour plot. We find that over the entire region exam¬ 
ined, G 15 remains large relative to kT, ranging from 
22 kT to 46 kT. G 15 monotonically increases with P along 
isotherms (Fig. fTT|) . and monotonically increases with T 
along isobars (Fig. fl5l). 

As shown in Fig. [TT] we note that the isotherms of 
G 15 are approximately linear in P at high P, but then 
start to curve as P passes through the vicinity of the Kt 
maximum. Indeed, there is reason to predict that these 
isotherms will pass through a minimum at negative P 
beyond the range of our current data. The density of the 
liquid and crystal are approaching one another as P de¬ 
creases, and can be expected to intersect in the negative 
pressure region. Near this intersection, the free energy 
cost to create ice nuclei of a given size should be lower 
than at higher or lower P, since under these conditions 
an ice embryo can form without a density fluctuation. 
This picture is also consistent with the requirement (im¬ 
posed by the Clapeyron equation) that the ice-liquid co¬ 
existence line must pass through a maximum T as P de¬ 
creases if the liquid and crystal densities intersect. Such a 
melting-line maximum has been observed for other water 
models; see e.g. Ref. [5Cl]. For a supercooled liquid near 
a melting-line maximum, the ice nucleation barrier must 


have a minimum as a function of P along an isotherm 
because the coexistence line (where the nucleation bar¬ 
rier diverges) is crossed by both increasing or decreasing 
P. 

Consistent with the above considerations, the onset of 
curvature in the isotherms of G 15 as P decreases through 
that of the Kt maximum reflects the change in liquid 
structure that occurs as the Widom line is crossed. As 
P decreases below that of the Kt maximum, the liquid 
enters the region where it becomes increasingly similar to 
the LDL phase, characterized by a high degree of local 
tetrahedral structure and a density approaching that of 
ice. Our data show that the onset of this crossover in the 
liquid phase is also reflected in the behaviour of the free 
energy of formation of small ice nuclei. 

The isobars of G 15 presented in Fig. [T5l show that G 15 
decreases with T at all P studied here. Fig. [15] also 
identifies the isobar corresponding to the estimated crit¬ 
ical pressure P c = 180 MPa of the LLPT for the ST2 
model, as well as the value of the critical temperature 
T c = 247 K EH- In the case that the supercooled liquid 
becomes unstable to ice formation, G* and n* should ap¬ 
proach zero, in which case G 15 should also approach zero. 
Our isobars of G 15 are too far away from T c to justify 
a quantified extrapolation. Nonetheless, our data sug¬ 
gests that future work to extend these curves to lower T 
would provide useful information. On the one hand, the 
slope of our isobar for P = P c suggests that a non-zero 
ice nucleation barrier exists in the vicinity of the criti¬ 
cal point of the LLPT, a result that is consistent with 
several recent free energy simulation studies of ST2 wa¬ 
ter [3, [50, HI, HH • On the other, our isobars for P < P c 
raise the possibility that the ice nucleation barrier may 
become low or even approach zero in the region of the 
phase diagram below the coexistence line of the LLPT, 
in the stability field of the LDL phase. If such a stability 
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limit for the LDL phase of ST2 water exists, clarifying its 
location would be valuable for understanding this phase 
of supercooled water, and for elucidating the influence of 
the LLPT on the ice nucleation process. 

V. DISCUSSION 

In summary, we have quantified the free energy of for¬ 
mation of small ice nuclei in the ST2 model of water. 
For the range of T and P studied here, we find that 
these small nuclei have the local structure of ice Ic, and 
that there is a significant free energy cost to form these 
ice nuclei in the vicinity of the Widom line. This result 
excludes the possibility that the liquid state of the ST2 
model is unstable to ice formation in this region of the 
phase diagram, and thus confirms that the thermody¬ 
namic anomalies previously identified in this region (e.g. 
the Kt maximum) occur in a well-defined metastable 
liquid state that is separated from the crystal by a sub¬ 
stantial free energy barrier. A similar conclusion was 
recently drawn for the case of the TIP4P/2005 model of 
water, though by different means pill|52|. 

From a methodological standpoint, our work illus¬ 
trates that two-dimensional umbrella sampling can be ex¬ 
ploited to quantify nucleation phenomena over a contin¬ 
uous range of T and P, and is a valuable approach when 
seeking to correlate thermodynamic behaviour and nu¬ 


cleation processes in simulations of supercooled liquids. 

For the specific case of supercooled water, our results 
show that we can detect the crossover from the HDL- 
like liquid to the LDL-like liquid by examining the nucle¬ 
ation behaviour of ice from the liquid phase. Since this 
crossover influences the free energy of formation of ice 
nuclei, we can also expect it to influence the nucleation 
rate of ice. This observation is significant because the 
principal challenge for the direct experimental character¬ 
ization of deeply supercooled water is the increasingly 
rapid conversion of the metastable liquid to ice as T de¬ 
creases. However, the ice nucleation process is strongly 
influenced by the nature of the metastable liquid in which 
the process begins, and our results demonstrate that ice 
nucleation can be used as a “probe” to reveal the proper¬ 
ties of the underlying supercooled liquid. Our work raises 
the possibility that an examination of ice nucleation rates 
in deeply supercooled water might be useful as a means 
to locate and study the thermodynamic anomalies of the 
supercooled liquid, including perhaps the LLPT itself. 
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